#!/bin/bash
#
# Copyright 2011 Google Inc. All Rights Reserved.
# Author: bst@google.com (Ben Stewart)
#
# Convert a BOM cyclone projection map to a format suitable for overlay on
# Google Maps.

SOURCE_URL='http://www.bom.gov.au/fwo/IDQ65001.gif'
TMP_BOM_IMAGE='qldcyclone.gif'
TMP_BOM_IMAGE_2='qldcyclone-edit.png'
TMP_GEOTIFF='qldcyclone.tiff'

# Output directory for the Maps tiles
TILE_DIR='tiles/'

# Background colour to remove.
COL_BG='#e6e6ff'
COL_PARALLEL='#70A270'

# Tile pyramid - min and max zoom levels
WND_ZOOM_MIN=3
WND_ZOOM_MAX=10

# Known points on the source image.
# These may need to be updated if the BOM changes their viewport of the image.
# The original image is assumed to be in Simple Cylindrical projection
# (EPSG:4326), with (X1,Y1) mapping to (LAT1,LNG1) etc. X,Y are in pixels with
# an origin at the top left of 0,0, and LAT,LNG are in degrees, with the origin
# at the prime meridian (Greenwich) and the equator.
X1=17
Y1=4
LAT1='-10.0'
LNG1='140.0'

X2=615
Y2=452
LAT2='-25.0'
LNG2='160.0'

# Extract a set of pixels in an image.
fetch_data () {
  FD_W="$1"
  FD_H="$2"
  FD_X="$3"
  FD_Y="$4"

  convert "${TMP_BOM_IMAGE}" -crop "${FD_W}x${FD_H}+${FD_X}+${FD_Y}" txt:/dev/stdout | \
    tail -n +2 | awk '// { print $3 }'
}

# Check for a crosshair of $COL_PARALLEL around ($1, $2)
bad_geom () {
  BG_X="$1"
  BG_Y="$2"

  COLDATA=$( (fetch_data 1 3 "${BG_X}" "$((BG_Y-1))";
      fetch_data 3 1 "$((BG_X-1))" "${BG_Y}") | uniq)

  test "${COL_PARALLEL}" != "${COLDATA}"
}

# Assert our settings are sane
if [ "$X1" -lt 0 ] || [ "$X2" -lt 0 ] || [ "$Y1" -lt 0 ] || [ "$Y2" -lt 0 ];then
  echo X and Y coordinates must be given in absolute pixel coordinates. >&2
  exit 1
fi

if [ "$WND_ZOOM_MIN" -lt 0 ] || [ "$WND_ZOOM_MAX" -lt 0 ] || \
  [ "$WND_ZOOM_MIN" -gt 20 ] || [ "$WND_ZOOM_MAX" -gt 20 ]; then
  echo Zoom levels must be in the range 0..20. >&2
  exit 2
fi

# Fetch the image from BOM
if ! curl -s -S -o "${TMP_BOM_IMAGE}" "${SOURCE_URL}"; then
  echo Could not fetch BOM source image. Giving up.
  exit 3
fi

# Verify the geometry
if bad_geom $X1 $Y1 || bad_geom $X2 $Y2 || bad_geom $X1 $Y2 || bad_geom $X2 $Y1; then
  echo It looks like the geometry specified is incorrect. The lat/long
  echo parallels are not appearing in the correct colour.
  echo ''
  echo Please check the downloaded imagery and update X[12], Y[12], LAT[12], LNG[12].
  exit 4
fi


# Convert the image to RGBA PNG with appropriate transparency so we can overlay
# on Maps.
if ! convert "${TMP_BOM_IMAGE}" -transparent "$COL_BG" \
  -fill 'transparent' -draw 'rectangle 416,433 627,446' \
  -draw 'rectangle 416,422 421,428' \
  -draw 'rectangle 513,422 531,428' \
  -draw 'rectangle 612,422 627,428' \
  -draw 'rectangle 500,452 542,458' \
  -alpha set -depth 8 -type TrueColorMatte -channel RGBA \
  -fill 'transparent' -opaque '#e8e1c4' \
  -fill 'transparent' -opaque '#70a270' \
  -fill 'transparent' -opaque '#9c988a' \
  -fill 'transparent' -opaque '#607960' \
  -fill 'transparent' -opaque '#9b9ba7' \
  -opaque '#c08000' \
  -fill 'rgba(192,192,192,0.4)' -opaque '#d3d3d3' \
  "${TMP_BOM_IMAGE_2}"; then
  echo Could not edit BOM image. >&2
  exit 5
fi

# Now that we have our transparent PNG, convert to GeoTIFF.
#
# TODO(bst): Big hack here. Assuming image field of view is constant. We
# really need a way to assert the pixels are correct in Imagemagick..
if ! gdal_translate -gcp "$X1" "$Y1" "$LNG1" "$LAT1" \
  -gcp "$X2" "$Y2" "$LNG2" "$LAT2" \
  -gcp "$X1" "$Y2" "$LNG1" "$LAT2" \
  -gcp "$X2" "$Y1" "$LNG2" "$LAT1" \
  -a_srs "EPSG:4326" \
  "${TMP_BOM_IMAGE_2}" "${TMP_GEOTIFF}"; then
  echo Could not translate map coordinates. Please check settings. >&2
  exit 6
fi

# Our GeoTIFF is ready for conversion from cylindrical project to Mercator.
if [ ! -d "${TILE_DIR}" ]; then mkdir "${TILE_DIR}"; fi
if [ ! -d "${TILE_DIR}" ]; then
  echo Could not create tile output directory. >&2;
  exit 7
fi

# Make sure tile data is empty - we don't want to serve stale tiles outside
# the current boundary.
echo Cleaning tile directory... >&2
find "${TILE_DIR}" -type f -delete || exit 8

# Invoke the tile generation process
echo Generating tiles. This may take 10-15 minutes, go make a coffee. >&2
gdal2tiles.py -z "$WND_ZOOM_MIN"-"$WND_ZOOM_MAX" "${TMP_GEOTIFF}" \
  "${TILE_DIR}" || exit 9

# Execute a mass compression for all the map tiles. This is required for
# serving the tiles.
echo Crushing PNG tiles. This will also take a while. >&2
find "${TILE_DIR}" -type f -name '*.png' -print0 | xargs -0 -P8 -I{} -L20 \
  pngcrush -e .png.out -brute -rem tEXt -rem tIME -rem iTXt -rem zTXt {} \
  || exit 10

# Blow away the big tiles, replace with crushed tiles.
find "${TILE_DIR}" -type f -name '*.png' -exec mv {}.out {} \; || exit 11

# TODO(bst): Extract the bounds for the overlay
# Transform to cylindrical projection so we can extract the image boundaries.
if ! gdalwarp -t_srs 'EPSG:4326' "${TMP_GEOTIFF}" \
  "${TMP_GEOTIFF%.tiff}-proj.tiff"; then

  echo Could not reproject the overlay to cylindrical projection to extract >&2
  echo the bounding box. >&2
  exit 12
fi

if ! (gdalinfo "${TMP_GEOTIFF%.tiff}-proj.tiff" | \
  awk '/^(Center)/ { lngC = $3; latC = $4;
           r = sub(/,$/,"",lngC); r = sub(/)$/,"",latC); };
       /^(Upper Left)/ { lng2 = $4; lat2 = $5;
           r = sub(/,$/,"",lng2); r = sub(/)$/,"",lat2); };
       /^(Lower Right)/ { lng1 = $4; lat1 = $5;
           r = sub(/,$/,"",lng1); r = sub(/)$/,"",lat1); };

       END { print "Please include the following in map_script.js:\n";
           print "initTileOverlay('\''LayerName'\'', '\''"              \
           "http://maps.gstatic.com/mapfiles/qldfloods/map-base-url/'\''" \
           ", '\'/\'', false, " lat1 ", " lng1 ", " lat2 ", " lng2        \
           ", '"${WND_ZOOM_MIN}"', '"${WND_ZOOM_MAX}"', ["                \
           latC ", " lngC ", '"${WND_ZOOM_MIN}"', '\''Place Name'\''"     \
           "], true);" } '); then

  echo Could not extract the bounding box for the overlay.
  exit 13
fi


